Synchronization and Collective Dynamics in A Carpet of Microfluidic Rotors 
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We study synchronization of an array of rotors on a substrate that are coupled by hydrodynamic 
interaction. Each rotor, which is modeled by an effective rigid body, is driven by an internal torque 
and exerts an active force on the surrounding fluid. The long- ranged nature of the hydrodynamic 
interaction between the rotors causes a rich pattern of dynamical behaviors including phase or- 
dering and self-proliferating spiral waves. Our results suggest strategies for designing controllable 
microfluidic mixers using the emergent behavior of hydrodynamically coupled active components. 

PACS numbers: 87.19.rh,07.10.Cm,47.61.Ne,87.80.Fe,87.85.Qr 



Introduction. Microorganisms and the mechanical 
components of the cell motility machinery such as cilia 
and fiagella operate in low Reynolds number conditions 
where hydrodynamics is dominated by viscous forces [1| . 
The medium thus induces a long-ranged hydrodynamic 
interaction between these active objects, which could lead 
to emergent many-body behaviors. Examples of such 
cooperative dynamical effects include sperms beating in 
harmony Q, metachronal waves in cilia formation 
of bound states between rotating microorganisms [6] , and 
flocking behavior of red blood cells moving in a capillary 
[7]. For a collection of free swimmers, such as microor- 
ganisms Q , hydrodynamic interactions have been shown 
to lead to instabilities |9|, llOll that can result in complex 
dynamical behaviors [l0|, In the context of simple 
microswimmer models where hydrodynamic interactions 
coupled to internal degrees of freedom can be studied 
with minimal complexity, it has been shown that the cou- 
pling could result in complex dynamical behaviors such 
as oscillatory bound states between two swimmers [l2| . 
and collective many-body swimming phases [HI, [HI . 

A particularly interesting aspect of such hydrodynamic 
coupling is the possibility of synchronization between dif- 
ferent objects with cyclic motions (3. [HI. fliMilj . This ef- 
fect has mostly been studied in simple systems such as 
two interacting objects or linear arrays and very little 
is known about possible many-body emergent behaviors 
of a large number of active objects with hydrodynamic 
coupling. For example, in a recent experiment [22|, Darn- 
ton et al. observed chaotic flow patterns with complex 
vortices above a carpet of bacteria with their heads at- 
tached to a substrate and their fiagella free to interact 
with the fluid (see also (lH). On the other hand, recent 
advances from micron-scale magnetically actuated tails 
[24j to synthetic molecular rotors [25| now allow fabrica- 
tion of arrays of active tails that can stir up the fluid. It 
is therefore very important to explore the possible com- 
plexity of the phase behavior of such an actively stirred 
microfluidic system. 

Here, we consider a simple generic model of rotors [26| 
positioned on a regular 2D array on a substrate and study 
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FIG. 1: (Color online). Schematic representation of the array 
of rotors. Inset: an immobilized bacterium with active fiag- 
ella as a possible realization of a rotor that can exert both a 
tangential drag and an active radial force on the fluid. 



their collective dynamics numerically. We find that the 
long-ranged hydrodynamic interactions could either en- 
hance or destroy ordering, depending on the degree of 
a built-in geometric frustration that originates from the 
interaction of the rotors with the fluid. More specifi- 
cally, our model adopts a fully synchronized state when 
the frustration is weak, and a randomly disordered state 
when it is maximally frustrated. Moreover, the dynam- 
ics of the system leads to self-proliferating spiral waves 
between the above two limiting behaviors. We also take 
into account thermal fluctuations of the rotors and map 
out the phase diagram of the system as a function of 
temperature and the degree of frustration. 

Model and Dynamical Equations. We consider an ar- 
ray of rotors that are assumed to be spherical beads of 
radius a moving on circular trajectories of radius 6, which 
are positioned on a rectangular lattice of base length d 
and at a height h above a substrate (see Fig. [1]). The i-th 
rotor is anchored at r^o to the surface of the substrate, 
which we take to be the x?/-plane. The instantaneous 
position of the rotating bead is = r^o + biii + he z , 
where the unit vector rii(t) = (cos fa (t) , sin fa (t) , 0) gives 
the orientation of the arm of the rotor. Because of 
the constraint that the bead is only allowed to move 
on the circular orbit of radius 6, the velocity of the 



rotor can be written as 
ti = e z x = (— sin fa, cos c 
gent to the trajectory. 



dt 



dt r * 



where 



i,0) is the unit vector tan- 



FIG. 2: (Color online). Snapshots of coarsening defects (with the greyscale representing cos0(r)) for (a) 6 = and spiral waves 
for (b) 5 = 45° and (c) 6 = 60°. These developed from random initial perturbations, (d) Spiral waves for 8 = 60° evolving 
from a defect pair, with a schematic picture of the director field (red, solid arrows) and the velocity field (blue, dotted arrows) 
near a +1 defect. 



We assume that the structure of the rotor is such that 
it drags the fluid with it as it moves (tangentially) along 
the circular trajectory, while it can also pump the fluid 
radially due to some internal degrees of freedom. The 
inset of Fig. [T] shows a possible realization of such a 
system in the case of bacteria whose heads are fixed on 
the substrate. In this example, the spinning rotation of 
the flagella would produce the pumping effect, while the 
precession of the axis of the flagella about the anchoring 
point would correspond to the tangential motion of the 
bead. Therefore, in our simplified model each rotor exerts 
a force, which can be decomposed into the radial, tangen- 
tial, and vertical components as = F n Ui-\-F t ti-\-F z e z . 
The velocity field of the fluid created by the rotors is 
given by v(r) = G(r — r$) • F^ where G(r) is the 
Blake-Oseen tensor [27[, which describes the hydrody- 
namic interaction near a flat surface with the non-slip 
boundary condition. Assuming that the arm length b 
and the height h are much smaller than the characteristic 
distance d between the rotors, we can use the 0(h 2 /d 2 ) 
approximation [4( , G a p(r) = f^^p- fc> r ot,f3 = x, y, and 
G az (r) = G za (r) = G zz (r) = 0, for a = x,y. Note that 
the z-component of the force is not coupled to flow and 
that the fluid velocity is lying in the x?/-plane. 

To obtain the flow velocity at the position of the rotors, 
we need to subtract the self-interaction, which involves 
the Stokes drag coefficient £ = GTrna. This yields 



dt 



37 ^2 — 



2tt 



(1) 



where r^- = 



Fj, ujt, n = Ft,n/((b) are the reduced 
forces, and 7 = Qh 2 jr\ = 6irah 2 is the hydrodynamic 
coupling constant. When the interaction is weak, we can 
simplify the phase equation (pp) following a standard pre- 
scription [28]. To this end, we rewrite it in terms of the 
slow variable <&i = fa — uj t t, and then integrate it over a 
cycle under the approximation that <$>i in the interaction 
term is constant over the period 2-k juj t . We obtain 



~dt I^T 2^ 
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■ sin(<l> i 



S), 



(2) 



uu f t uJ 2 . This equa- 
% [3^ . In this form, 



where S = tan -1 (co>t/co> n ) and uj = \/uj+ 
tion is correct to 0(j/(d 3 sin 5)) [2£ 
the hydrodynamic coupling becomes isotropic and our 
system resembles existing models of non-locally coupled 
oscillators with phase delay 29|, 30| . 

Simulation Method. The model is implemented on a 
L x L square lattice with the grid size and the phase 
equation Eq. (pQ) is solved by Euler method with the 
time step At. The system size used is L = 256 for most 
of the results below, while we have also used L = 128 for 
obtaining some of the statistical data. 

We have imposed periodic boundary condition and 
computed the velocity field at every time step by Fourier 
transformation. For 5 > 0, we also solved the reduced 
phase equation Eq. (j2j), to compare with the solution 
of Eq. (pp), and found a very good agreement. We have 
also incorporated thermal fluctuations, by adding an un- 
corrected Gaussian noise Qi(t) to the RHS of Eq. (pp). 
The noise is assumed to have zero mean and its fluc- 
tuations are controlled by the rotational diffusion con- 
stant of the bead D r = k B T/((b 2 ) as (^(^(t 7 )) = 
2D r Sij5(t — t'). We define the reduced (effective) tem- 



perature r 



D r d 3 



k B Td 3 



-. For typical values 



juj 367r 2 r]a 2 b 2 h 2 uj 
oia^b^h^l //m, d ~ 10 //m, uj ~ 10 2 Hz, with 
r] = 1 x 10" 3 Pa s and k B T = 4 x 10" 21 J, we have 
r ~ 10 _1 . Note that the sharp dependence of r on a, 
6, h, and d makes it easy to control the reduced temper- 
ature by changing the size or density of rotors. In our 
simulations, we have used the parameter values 7 = 0.1 
and uj = 0.1, with d = 1 and At = 0.1. First we turn off 
the thermal noise (r = 0) and vary the force angle 6 to 
study the pattern dynamics. 

Pumping- driven Rotors. When 5 = 0, each rotor 
pumps the fluid radially, and is driven by the fluid flow 
generated by the other rotors. In this case, the initial ran- 
dom perturbations develop into topological defects (sin- 
gularities of the phase field far) of winding numbers ±1, 
which coarsen by collision of +1- and —1- defects and fi- 
nally disappear to establish global synchronization. Fig- 
ure [2a) shows a snapshot of the coarsening defects at 
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FIG. 3: (Color online), (a) Correlation length £ as a function of time and for different values of 8. For 5 — 0° it is well fitted 
by £ oc t 0-75 . For S < 5 C — 40° the correlation length diverges after t = 10000 (not shown), while for S > 5 C it remains finite, 
(b) Equilibrium order parameter S as a function of temperature r and for different values of 6. (c) Phase diagram and contour 
of the correlation length £. The correlation length is smaller for larger 6 and higher temperature. 



t = 10000. To characterize the phase ordering dynam- 

1/2 

ics, we define the correlation length £ = ((V0(r)) 2 ) r , 
where V0 is the spatial gradient and (...) r means spatial 
average. As shown in Fig. [3fa), we find that £ as a func- 
tion of time is well fitted by the power law £ ex t v , with 
v = 0.75. The scaling of the hydrodynamic interaction 
G(r) ~ | r | 3 and Eq. (pQ) suggest that the characteristic 
timescale of a pattern is proportional to its size, which 
would mean v = 1 . The difference between the numerical 
and scaling exponents suggests violation of dynamic scal- 
ing, which is characteristic to coarsening of point defects 
in two dimensions (3l| . 

Torque-driven Rotors. When S = 90°, each rotor is 
driven by an active torque and exerts a force tangential 
to its orbit. In this case, we find that the system reaches 
a disordered state in which spatial correlation is almost 
completely lost. The absence of orient ational correlation 
can be understood as follows. The average flow created 
by a rotor is perpendicular to its arm, and a neighbor- 
ing rotor tends to align with the flow. Thus the two 
rotors tend to be perpendicular to each other on aver- 
age. However, it is not possible that every pair of rotors 
have their arms vertically crossed (geometric frustration), 
and hence the system evolves towards randomly oriented 
states. 

Rotors driven by Pumping and Torque. In the gen- 
eral case, the rotor is driven by an active torque while it 
pumps the fluid radially, and the total force exerted on 
the fluid has an angle 0° < 6 < 90° with respect to the 
arm of the rotor. We varied the parameter 5 and found 
two types of dynamical behavior. For 0° < S < 40°, 
the globally synchronized state is still obtained as the 
final state. However, for 40° < S < 90°, we find that 
the dynamical steady state of the system involves self- 
proliferating spiral waves as shown in Fig. EJb) and (c). 
Moreover, we find that the correlation length con- 
verges to a finite value as t —> oo as shown in Fig. EJa), 
and that the equilibrium correlation length decreases as 
S is increased. 



The flow pattern is locally correlated with the rotor's 
director n(r) = (cos 0(r), sin0(r)). The surface flow ve- 
locity v(r) makes the angle S with n(r) except at the 
core of the defect (see the movies [32] for comparison of 
the two fields). This observation leads us to an intuitive 
interpretation of the spiral waves. In the vicinity of a 
+1 defect from which the director emanates radially, the 
rotors create an outgoing flow that has an anti-clockwise 
slant with respect to the radial direction, and hence form 
an anti-clockwise spiral; see the inset of Fig|2jd). The 
spiral is tighter for a larger force angle 5. We confirmed 
this scenario by choosing a defect pair as the initial con- 
figuration and following its evolution; see Fig|2fd) and 
the corresponding movie [32]. Initially, clockwise and 
ant i- clockwise spirals are formed around —1- and +1- 
defects, respectively. Then the director is randomized on 
the thinning spiral arm, which collapses and proliferates 
a cascade of new defects. 

Thermal Fluctuations. We then introduce the ther- 
mal torque and study the phase behavior of the system. 
In Fig. [3f b), we plot the equilibrium order parameter 
S = | (cos 0)| as a function of the effective temperature. 
For 6 = 0, we find a critical temperature r c at which S 
vanishes as S ex yV c — r to a good approximation (33| . 
We find the critical temperature for this phase transi- 
tion as r c = 0.76, which is about 30% smaller than the 
mean- field value r c = 1.08 by Guirao and Joanny [4]. 

As we increase the phase delay S up to S c = 40°, the 
critical temperature is lowered down to r c = 0.61. For 
S > 5 c, the order parameter S is less than 1 even at 
r = and is smaller for a larger system size L, suggest- 
ing that S = in the thermodynamic limit. However, 
the existence of local order (spiral waves) is reflected in 
the finite- L data, using which we can define the transition 
temperature r c in the same way as mentioned before [33[ . 
The resulting phase diagram is shown in Fig. [3fc). We 
distinguish three regions: (O) ordered, which is a distinct 
thermodynamic phase characterized by global synchro- 
nization, and (S) spiral waves and (D) disordered, which 
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are inherently the same phase but with different local or- 
dering and dynamical structure. Also shown in Fig. [3fc) 
are the contours of the correlation length £. While the 
O-D and O-S transitions are sharp, the S-D transition is 
a crossover characterized by gradual decrease of £. We 
also note that the frustrated state for S = 90°, r = 
and the thermally-disordered state for S — 0°, r > r c are 
qualitatively different, though they are not distinguished 
in the phase diagram. 

Discussion. The case of no active torque (S = 0) 
could be regarded as a simplified and idealized model 
of bacterial carpets [lH, [23| . Our model reproduced the 
enhancement of orientational ordering, while it predicts 
global ordering and not the finite-size correlation as ob- 
served in the experiments. The experimental patterns 
might be explained by some kind of frozen disorder in 
the flagellar configuration, which can be readily incor- 
porated in our model 0, EHj. The case of S = 90° is 
realized by rigid spheres without pumping. It is related 
to a recently studied model of two rigid spheres making 
tilted elliptic orbits [17j, which show both in-phase and 
anti-phase synchronization. Our results suggest that the 
interaction between many of such rotors is frustrated and 
the system does not attain full synchronization. Spiral 
waves for finite phase delay S have been observed in pre- 
vious models of 2D coupled oscillators {2j| [lo|, [36j . How- 
ever, in our case, the pattern is intrinsically turbulent 
and self-proliferating, in contrast to the case of finite- 
range coupling, for which t he p hase is spatially smooth 
except near the defect core 

In conclusion, we have introduced a generic model of 
microfluidic rotors that shows a variety of dynamical pat- 
terns including global synchronization, fully disordered 
states, and self-proliferating spiral waves. The patterns 
are sensitively controlled by the angle of active force 5 
(the degree of frustration) and the temperature r. Our 
results suggest that arrays of active microfluidic compo- 
nents could be designed to induce a rich variety of dy- 
namical behaviors in the vicinal fluid, and could be used 
to make switchable microfluidic mixers. 
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